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ABSTRACT 


A numerical approach is presented for design sensitivity analysis. 

The approach is based on perturbing the design variables and then using 
iterative schemes to obtain the response of the perturbed structure. A 
forward difference formula then yields the approximate sensitivity. 
Algorithms for displacement and stress sensitivity as well for eigenvalue 
and eigenvector sensitivity are developed. Results for the stress 
sensitivity problem are compared with the semi-analytical method. Examples 
are considered in structures and fluids. 

INTRODUCTION 

Iterative methods are presented for obtaining design sensitivity 
coefficients (or derivatives) of implicit functions. Design derivatives 
are important not only in gradient-based optimization codes, but also for 
examining trade-off's, system identification, and probabilistic design. 
Iterative methods are presented for both the algebraic and eigenvalue 
problems; stress, eigenvalue and eigenvector derivatives are considered. 

The iterative approaches provide approximate derivatives. They are very 
simple to implement in a program, especially for calculation of eigenvector 
derivatives. The idea of using iterative methods for a class of problems 
was suggested for one dimensional problems in (ref. 1). Here, this idea is 
developed to handle the matrix algebraic as well as the generalized 
eigenvalue problems. 

The basic idea behind the approach is as follows. Let 

9 = g(b,£) (1) 

be a continuously differentiable function of a design variable vector b of 
dimension (kxl), and a state variable vector y of dimension (nxl). The 
state variables are implicitly dependent on design through the n state 
equations of the form 


4(b ,£) - 0 (2) 

Let b° be the current design and y° be the associated state variable 
vector. The problem of concern is to find the sensitivity, dg/db, at the 
current design. The iterative method is based on perturbing each" design 
variable, in turn, as 


b- - b° + e (3) 

Equation (2) now becomes 

4(b £ ,Z £ ) - 0 (4) 

Now, a modified residual-correction or Newton-Raphson technique is 

applied to solve Eq. (4), treating y e as the vector of unknowns. 

Then, the sensitivity vector is given approximately by 

dg/db. - [g(b E , /)-g(b°,y°)]/e (5) 
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For the eigenvalue problem, as discussed later, the system in Eq. (4) 
is augmented by a certain orthogonality relation. Note that certain 
coefficient matrices involving stiffness, mass, etc. have been decomposed 
at the current design while solving for ^°. The iterative approach 
presented here can be viewed as re-analysis schemes used to solve Eq. (4), 
which uses the already decomposed matrices. Since the perturbation e is 
very small, the iterative schemes converge very rapidly. 

DISPLACEMENT AND STRESS SENSITIVITY 

A finite element model of the structure is assumed. The problem of 
obtaining design derivatives of displacements and stresses Is now 
considered. Consider a function 


9 - fl(b,z) (6) 

which represents a stress constraint, with b - (kxl) design vector and z - 
(nxl) displacement vector which is obtained~from the finite element 
equations 

K(b) z - F(b) (7) 

where K is an (nxn) structural stiffness matrix, and F is an (nxl) nodal 
load vector. Let b° be the current design. At this ?tage, the analysis 
has been completed. Thus, the decomposed K(b)° and z° are known. 

The derivative of the function g with respect to the ith design 
variable is given by 


dg/db-j - dg/dbi + dg/dz • dz/db^ (8) 

The partial derivatives dg/db and dg / dz are readily available using the 
finite element relations. The problem, therefore, is to compute the 
displacement sensitivity, dz/db. An iterative approach for computing this 
is now given. 

Corresponding to the ith design variable, let the perturbed design 
vector, b e , be defined as 

b e - (b{, b° b°. + e b°) T (9) 

The perturbation e is relatively small, and a value of 1% of b-j has found 
to work well in practice. The choice of e is based on balancing the 

truncation and cancellation errors. The problem is to find z e , the 
solution of 

K(b e ) z e = F(b e ) (10) 

using the decomposed K(b°) and z°. A modified version of the 
residual-correction scheme given in (ref. 2) is given below. 
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Algorithm 1 (Displacement and Stress Sensitivity) 


Step (0). 

Set j-0. Choose the perturbation e 
A. 

and the error tolerance 


Define b e as in Eq. (9). 


Step (i). 

Calculate the residual rJ from 


Step (ii). 

- K(b e ) z J - F(b e ) 
Solve for the correction eJ from 

(11) 


K(b°) eJ - -rJ 

(12) 

Step (iii). 

Update 



zJ + l - zJ + eJ 

(13) 

Step ( i v ) . 

Check the convergence criterion 



1 1 zJ + l - zJ|| s A 

(14) 


If (14) is satisfied, then set z e - zJ +1 and compute the 
displacement sensitivity as 


dz / db-j - (z e - z°) / e (15) 

The stress sensitivity can be recovered from Eq. (8). If Eq. 
(14) is not satisfied, set j * j+1 and re-execute steps 
(i)-(iv) above. 

Numerical results and comparison with the exact and semi-analytical 
methods discussed in the literature are presented subsequently. 
Theoretically, it can be shown that the above scheme will converge provided 
[ 2 ]: 


r 0 Cl - K'V)K(b e )] < 1 (16) 

where r 0 ( A) =* spectral radius of the matrix A, which is the maximum size 

of the eigenvalues of A. In the problem considered here, K(b°) and K(b e ) 
are roughly equal owing to e being small, and (16) can generally be 
expected to hold. 


EIGENVALUE AND EIGENVECTOR SENSITIVITY 

Eigenvalue sensitivity is useful when resonant frequencies or critical 
buckling loads need to be restricted. Exact analytical expressions for 
eigenvalue sensitivity can be readily derived for the case of non-repeated 
roots [3]. The problem of obtaining eigenvector sensitivity, on the other 
hand, is more complicated and is an area of current interest [4-7]. 
Eigenvector sensitivity is useful in obtaining the design derivatives of 
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forced dynamic response. Here, an iterative approach is presented for 
approximate derivatives of eigenvalues and eigenvectors. The approach is 
particularly easy to implement in a program and provides both eigenvalue 
and eigenvector derivatives simultaneously. Further, the derivative of a 
particular eigenvector does not require knowledge of all eigenvectors of 
the problem, as with certain analytical methods. 

Consider the generalized eigenvalue problem 

K(b)y. - A M(b )i (17) 

where A is a particular non-repeated eigenvalue and is the associated 
eigenvector. For the frequency problem, K and M represent the structural 
stiffness and mass matrices, respectively. For the buckling problem, K and 
M represent the structural stiffness and geometric stiffness matrices, 
respectively. It is desired to find the sensitivities dA/db and dy/db. 

Let b° be the current design vector and (A 0 , y°) be a given 

eigenvalue-eigenvector pair at the current design. Let b e be a perturbed 
design vector as given in Eq. (9). The residual is given by 

R - K(b e )y e -A £ M(b e )y e (18) 

The object is to solve the nonlinear equations R - 0 for the unknowns A £ and 
y e ; the Newton-Raphson technique is used for this purpose. The Jacobian J 

of the system in Eq. (18) is [dR/dy e , 3R/5A ]. The Newton-Raphson equations 
are consequently: 


[K(b e H £ M(b e ) -M(b e )y ] = -R (19) 

Note, however, that Eq. (19) represents a system with n equations and (n+1) 
unknowns; an additional equation is needed. This additional equation is 
obtained by introducing the normalization condition 

-£ t M 6y. = 0 (20) 

which states that the change in the eigenvector is orthogonal to the 
original eigenvector with respect to the mass matrix. In fact, the above 
scheme has been used as a re-analysis approach in (ref. 8). Here, an 
additional modification is made: the Jacobian matrix in Eq. (19) is 

modified by replacing K(b e ) by K(b°), y e by y° and A £ by A 0 . The 

motivation for this, as in the previous section, is tc preserve a constant 
coefficient matrix in the iterative scheme. The resulting efficiency has 
not been found to affect the convergence of the procedure owing to the 
relatively small size of e. The above modifications lead to an iterative 
scheme based on solving the system. 

® & ■ (i) 
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where 


K(b°)-A 0 M(b°) 


-y° M(b°) 



( 22 ) 


The coefficient matrix is symmetric and nonsingular for the case of 
non-repeated roots [8]. Gaussian elimination can be used to solve Eq. (21). 
The algorithm for eigenvalue-eigenvector sensitivity is now given. 


Algorithm 2 (Eigenvalue-Eigenvector Sensitivity) 

Step (0). Set j-0. Choose the perturbation e and the error tolerances 
Aj and A£. 

Define b e as in Eq. (9). Decompose the matrix C given in Eq. 

( 22 ). 

Step (i). Define the residual 


Step (ii). 


Step ( i i i ) . 


Step (iv). 


RJ = K(b e )yj - A jM(b e )yJ 


Solve the algebraic equations 



for <5y and <5A. 

Update 

yj+1 - yj + <5y 
Aj+i ■ Aj + 6A 
Check the convergence criterion 


6y 1 1 i Aj , | 6A| i A 2 


(23) 


(24) 


(25) 

(26) 


If (26) is satisfied, then set y e =yJ + l, A £ = Aj+i and 
compute the sensitivity as 


dy/db-j « (y e -y°)/e 
dA/db-j = (A e -A 0 )/e 


(27) 
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If (26) is not satisfied, set j = j+1 and re-execute steps 
(i)-(iv) above. 

Numerical results are presented in the following section. 

NUMERICAL RESULTS 


Thin plate problem 

Consider the plane stress problem in Fig. 1, where inverse thicknesses 
are the design variables. That is, the reciprocal of the plate thickness 
is chosen as a design variable. Inverse design variables are used in 
optimal design literature because they linearize the stress function and 
lead to improved convergence. The stress constraint function is the 
von Mises failure criterion in element j, given by 

9j - ovM/®a _1 (28) 

where ovm 2= o x 2 + °y 2 "* °x°y + 3x X y2 and o a - constant allowable stress 
limit. Constant strain triangular elements are used. For brevity, only 
the design sensitivity coefficients, dg^g/dbig and dg 24 /db 24 , are presented 
in Table 1. The sensitivity vectors have been obtained using Algorithm 1 
discussed earlier. In Table 1, the results obtained by the iterative 
method are compared with the semi-analytical method used widely in the 
literature, based on the formula in Eq. (8) with dz/db-j obtained from 

K(b e )-K(b°) o [f(b e )-F(b 0 )] 

K dz/db i - - - z + (29) 

The results are also compared with the exact sensitivity obtained using 
analytical derivatives. It is interesting to note that the semi-analytical 
method yields the same result as the first iteration of the iterative 
method. However, the iterative method further improves upon this and 
approaches the exact sensitivity (Table 1). While all methods yield values 
of acceptable accuracy, the comparison serves to illustrate the nature of 
the iterative process. This aspect is shown graphically in Fig. 2. It is 
noted that when using direct variables (as opposed to reciprocal 
variables), the semi-analytical method yields essentially exact sensitivity 
owing to the fact that stiffness is a linear function of design variables. 



L i ( Poisson ratio ) * 0.2 

o ( Allowable stress ) - 20000. psi 


Figure 1. 
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Exact 



Iterative 5 
Iterative 4 
Iterative 3 
Iterative 2 
Iterative 1 
= Semi-analytical 


» — »,• 

i 

b°+e 

1 

Figure 2. 
Table 1. 

i 

Method 

dg 19 / db xg 

dg24 /db 24 

1 

8.7098 

5.6437 

Iterative 2 

8.7949 

5.6980 

3 

8.7957 

5.6986 

Semi-analyti cal 

8.7098 

5.6437 

Exact 

8.7969 

5.7002 


Plane Frame Problem 


Consider the frame structure in Fig. 3. The design variables 
associated with the I-section are b = (h, w, t w , tf) as shown in Fig. 3. 

The current design is b = (3.0, 3.0, 0.3, 0.5) for each element. The 
sensitivity of the lowest eigenvalues and corresponding eigenvector 
obtained using Algorithm 2 given earlier is presented in Tables 2 and 3, 
respectively. For the eigenvector, only selected sensitivity coefficients 
are presented for brevity. The maximum number of iterations required for 
an error tolerance of 10"? is five. Thus, we see that convergence of the 
algorithm is very rapid and simple to implement. Also, the algorithm does 
not require computation of all eigenvectors to find the sensitivity of a 
few specific eigenvectors. However, if the sensitivity of all eigenvectors 
is required, then alternative approaches may be preferable. 


720 











Table 3. 


No. of degree of freedom 

dy/dbs 

dy/dbi7 

4 

0.036758 

0.021438 

5 

0.000370 

0.002413 

6 

-0.000868 

0.003890 

7 

-0.028564 

-0.007628 

8 

0.000378 

0.002547 

9 

-0.002667 

-0.000389 

10 

-0.028327 

-0.007628 

11 

-0.000653 

-0.002547 

12 

0.001614 

-0.000389 

13 

0.036339 

0.021438 

14 

-0.000407 

-0.002413 

15 

0.000147 

0.003890 


Fluid Mechanics Problem 


The objective of this problem in Fig. 4 is to obtain the sensitivities 
of the maximum absolute eigenvalue and eigenvector of the amplification 
matrix G of the incompressible Euler equations in fluid mechanics (ref. 9). 
This problem is motivated from a study of the stability of the 
computational algorithm. The Euler equations are 


[I - At D + -1 (l-cos0 x ) 1+1 jg A sin0 x ] (I - At D)' 1 

+[I - At D + -j (l-cos0 y ) 1+1 0 B sin0 y ] (I - At Of 1 

+[I - At D + ^ (l-cos0 z ) 1+1 ^ C sin9 z ] (G - I) 

- At D [(1 - cos0 x ) 2 + (1 - cos0 y ) 2 + (1 - cos0 z ) 2 ] I 

- 1 Ifth sin0 x + sin0 y + sin9 z^ 

where, l is a (4x4) identity matrix. 

The source Jacobian Matrix is 
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Data 


Grid sizes in x, y and z directions are Ax = 1/16, Ay - tr/32, Az * 1/32. 
Parameter of time-derivative term is 6 * 1. 

Radius is r-2, angular velocity of propeller is oi-2. 

Parameter CFL (Courant-Friedrichs-Lewy Number) is - 5. 

Fluid velocities in x, y, z directions are u x » 0.5, Uy - 1, u z - 1. 

The implicit second-order artificial viscosity is e-j - 0. 

The explicit fourth-order artificial viscosity is e e - 0. 

The lower boundaries of wavenumbers 0 X , 0y, 0 Z are = 0, and the upper 
boundaries are = n. 

Results 


The optimum values of wavenumbers 0 X , 0y, 0 Z at which the absolute 
value of maximum eigenvalue are maximum, are obtained by using the 
optimization program LINRM [10]. The results are 0 X - 0y * 0 Z - ir/2. All 
sensitivity calculations are now done at these values or 0 X , 0y and 0 Z . 

The sensitivity of the maximum eigenvalue, absolute maximum eigenvalue and 
corresponding eigenvector are shown in Tables 4 and 5. 



z 



Figure 4. 
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Table 4. Eigenvalue Sensitivity for Fluid Mechanics 
Problem 


A - 1 x 10" 8 , e - 0.01 


Design Variable 

Eigenvalue Sens 

itivity 

dA/db 

d| A | /db 


-0.03928 -0.11272 i 

-0.08543 

Uy 

0.25361 -0.02037 i 

0.21790 

“z 

-0.23827 -0.36236 i 

-0.37489 

CFL 

0.13254 +0.06004 i 

0.14541 

£i 

-0.12951 -0.29317 i 

-0.24668 

£e 

-0.84451 -0.46701 i 

-0.54753 

to 

-0.01807 +0.02823 i 

-0.00358 

e x 

-0.00138 -0.00002 i 

-0.00125 

e y 

-0.00182 -0.00010 i 

-0.00167 

e z 

-0.00204 -0.00189 i 

-0.00267 


bo =. (0.5, 1, 1, 5, 0, 0, 2, n/2, n/2, ir/2) 
A° = 0.94882 + 0.47287 i 
| A 0 | - 1.0601 


Table 5. Eigenvector Sensitivity for Fluid Mechanics Problem 


A - 1 x 10"8, e = 0.01 

I '..i f ■ 


Degree of Freedom 

dy/d CFL 

dy/du 

1 

-0.22740 -0.27610 i 

0.03483 +0.01911 i 

2 

-0.34796 -0.03808 i 

0.00969 -0.00969 i 

3 

-0.20260 -0.10825 i 

0.07636 +0.00170 i 

4 

-0.19075 -0.10196 i 

0.03520 -0.01693 i 



SUMMARY AND CONCLUSIONS 


A numerical method has been presented for design sensitivity analysis. 
The idea is based on using iterative methods for re-analysis of the 
structure due to a small perturbation in the design variable. A forward 
difference scheme then yields the approximate sensitivity. Algorithms for 
displacement and stress sensitivity as well as for eigenvalues and 
eigenvector sensitivity are developed. The iterative schemes have been 
modified so that the coefficient matrices are constant and hence decomposed 
only once. The convergence is found to be very rapid. Further, 
implementation of the algorithms is simple. 
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